#!/usr/bin/env python
"""

    nth_percentiles.py
    [--log_file PATH]
    [--verbose]

"""

################################################################################
#
#   nth_percentiles
#
#
#   Copyright (c) 7/8/2010 Leo Goodstadt
#
#   Permission is hereby granted, free of charge, to any person obtaining a copy
#   of this software and associated documentation files (the "Software"), to deal
#   in the Software without restriction, including without limitation the rights
#   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#   copies of the Software, and to permit persons to whom the Software is
#   furnished to do so, subject to the following conditions:
#
#   The above copyright notice and this permission notice shall be included in
#   all copies or substantial portions of the Software.
#
#   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
#   THE SOFTWARE.
#################################################################################

import sys, os


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Functions


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#_________________________________________________________________________________________

#   nth_percentiles

#_________________________________________________________________________________________
ROUNDING      = 0
AVERAGING     = 1
NEAREST_EVEN  = 2
KNOTS_INTERPOLATE = 5
S_INTERPOLATE = 7
def nth_percentiles (array, N = 5, method = S_INTERPOLATE):
    """
    Algorithm From R quantile (Discontinuous sample quantile)
        Histogram of underlying distribution quantiles
        Weighted averages of consecutive order statistics

    method = ROUNDING:
        Rounding down to find nearest order statistic
    method = AVERAGING:
        Averaging at discontinuities
    method = AVERAGING:
        Nearest even order statistic (SAS definition)
    method = KNOTS_INTERPOLATE:
        Interpolation where middle values (excluding 0% / 100%) are values midway
        through the steps of the empirical cdf.
    method = S_INTERPOLATE:
        Interpolation used by S
        Evenly spaced sampling of cdf including 0% / 100%
        Sampling intervals of the cdfs are therefore 100/N% for N-tiles
    """
    if len(array) == 0:
        return [0] * (N + 1)
    s_array = sorted(array)
    sz = len(array)
    last_pos = sz - 1

    # first value in ordered statistic
    reduced_values = [s_array[0]]


    if method == ROUNDING:
        for i in range(1, N):
            np = i * 1.0 * sz / N - 1.0
            pos = int(np)
            if np - pos != 0:
                pos += 1
            reduced_values.append(s_array[pos])
    elif method == AVERAGING:
        for i in range(1, N):
            np = i * 1.0 * sz / N - 1.0
            pos = int(np)
            if np - pos == 0:
                reduced_values.append((s_array[pos] + s_array[pos + 1]) * 0.5)
            else:
                reduced_values.append(s_array[pos + 1])
    elif method == NEAREST_EVEN:
        for i in range(1, N):
            np = i * 1.0 * sz / N - 0.5
            pos = int(np)
            gamma = np - pos
            if gamma == 0 and pos % 2 ==0:
                pos = pos - 1
            reduced_values.append(s_array[pos])
    elif method == S_INTERPOLATE:
        for i in range(1, N):
            p = i * 1.0 / N
            np = sz * p
            pos_frac = np + 1 - p
            pos      = int(pos_frac)
            gamma    = pos_frac - pos
            val = (1 - gamma) * s_array[pos - 1] + gamma * s_array[pos]
            reduced_values.append(val)
    elif method == KNOTS_INTERPOLATE:
        for i in range(1, N):
            np = sz * i * 1.0 / N  + 0.5
            pos      = int(np)
            gamma    = np - pos
            val = (1 - gamma) * s_array[pos - 1] + gamma * s_array[pos]
            reduced_values.append(val)


    else:
        raise Exception("Unknown method %s" % str(method))

    # last value in ordered statistic
    reduced_values.append(s_array[-1])
    return reduced_values


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Main logic


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
if __name__ == '__main__':
#   debug code not run if called as a module
#
    #88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

    #   Testing


    #88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

    import unittest
    from random import randint, shuffle
    class Test_nth_percentiles(unittest.TestCase):

        #       self.assertEqual(self.seq, range(10))
        #       self.assert_(element in self.seq)
        #       self.assertRaises(ValueError, random.sample, self.seq, 20)
        #
        @staticmethod
        def create_distribution (cnt_values, max_value):
            import random

            # sort into clusters
            values = []
            while len(values) < cnt_values:
                test_value = max_value
                while test_value >= max_value:
                    test_value = int(max_value * random.lognormvariate(1, 5) * 0.5)
                values.append(test_value)

            #values = range(10)
            print "  cnt values = ", len(values)
            print "  max value  = ", max(values)
            return values

        def test_nth_percentiles(self):
            """
                test pairwise tuples
            """
            values = self.create_distribution(100, 1000)
            n_quantiles = 8
            quartile_sz = 100.0 / n_quantiles
            for (i, (v1, v2, v3, v5, v7)) in enumerate(zip( nth_percentiles (values, n_quantiles),
                                                nth_percentiles (values, n_quantiles, AVERAGING),
                                                nth_percentiles (values, n_quantiles, NEAREST_EVEN),
                                                nth_percentiles (values, n_quantiles, KNOTS_INTERPOLATE),
                                                nth_percentiles (values, n_quantiles, S_INTERPOLATE),
                                                )):
                print "  %3d%%ile" % (i * quartile_sz), "%6s" % v1, "%6s" % v2, "%6s" % v3, "%6s" % v5, "%6s" % v7
            #print nth_percentiles (values, 10)
            #print nth_percentiles (values, 20)
            #print nth_percentiles (values, 100)


    #
    #   call unit test without parameters
    #

    if sys.argv.count("--debug"):
        sys.argv.remove("--debug")
    unittest.main()







